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For general two-electron two-centre integrals over Slater-type orbitals (STOs), the 
use of the Neumann expansion for the Coulomb interaction potential yields infinite 
series in terms of few basic functions. In many important cases the number of terms 
necessary to achieve convergence by a straightforward summation is large and one 
is forced to calculate the basic integrals of high order. We present a systematic ap¬ 
proach to calculation of the higher-order terms in the Neumann series by large-order 
expansions of the basic integrals. The final expressions are shown to be transpar¬ 
ent and straightforward to implement, and all auxiliary quantities can be calculated 
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appropriate convergence accelerators. 

Keywords: Slater-type orbitals, integrals, large-order expansion, Neumann series 


^Electronic mail: lesiuk@tiger.chem.uw.edu.pl. 


1 





I. INTRODUCTION 


From a purely theoretical point of view, Slater-type orbitals^ (STOs) constitute a more 
convenient basis set for calculations in molecular physics than the widely used Gaussian-type 
orbitals^ (GTOs). In fact, unlike GTOs, STOs are able to satisfy the Kato’s cusp condition^ 
at the electron-nucleus coalescence points and their exponential decay at large electron- 
nucleus distances coincides with the asymptotic form of the electronic density^ (if nonlinear 
parameters are suitably chosen). Only severe difficulties in calculation of the electron repul¬ 
sion integrals made the use of STOs drastically limited. Nonetheless, a considerable interest 
remained in this field^ - —. 

In recent three paper series^— calculation of the STOs integrals has been reconsidered 
and new analytical or seminumerical methods for their computation have been proposed. 
This allowed to perform calculations for the beryllium dimer with STOs basis sets up to 
sextuple £ quality, reaching the so-called spectroscopic accuracy (few wavenumbers, cm -1 ). 
Additionally, it was found that the Coulomb, (aa\bb), and hybrid, (aa|a&), integrals are not 
troublesome and are computed with a decent accuracy and speed for a reasonable range 
of nonlinear parameters (and quantum numbers). Calculation of the exchange integrals, 
(ab\ab), is more involved. The Neumann expansion of the interaction potential, which is 
the method of choice, gives rise to infinite series. In many important cases the required 
accuracy is obtained after summing 20-30 terms. However, there are situations where a larger 
number of terms is necessary to achieve convergence which makes calculations significantly 
more expensive. This is one of the major reasons for the STOs vs. GTOs gap in the 
computational timings. 

It seems reasonable to expect that the higher-order terms in the Neumann expansion 
do not need to be computed with general techniques but a suitable large-order expansion 
can be devised. This would allow to reduce the computational burden significantly, as the 
asymptotic expansions of such kind are typically more robust than the general expressions. 
Therefore, the main purpose of this paper is to derive systematic large-order approximations 
of all basic quantities appearing in the Neumann expansion of the STOs exchange integrals, 
and provide necessary numerical tests. Resulting expressions can be readily incorporated 
into existing STOs integral codes. 

Since the present paper is concentrated solely on the Neumann expansion of the inter- 
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action potential with application to the STOs electron repulsion integrals, a brief survey of 
the literature on this topic is mandatory. Relevant mathematical details will be given in 
the next section. Possibly the first method utilising the Neumann expansion was reported 
by Ruedenberg^i^ who introduced general expressions based on charge distributions of 
both electrons. Later, this approach was extended by applying a straightforward numerical 
integration^ 1 ^ to avoid several difficult analytic rearrangements. Kotani^ provided many 
tools and expressions enabling fully analytical (albeit recursive) techniques to be used. Re¬ 
cursive approach was later pursued by Harris^ who invoked the theory of spherical Bessel 
functions to simplify the existing theory and discovered many useful additional relations. A 
considerable interest remained in the field despite GTOs were clearly taking over the role 
of routine basis set in quantum chemistry. Many changes were introduced in how individ¬ 
ual terms in the Neumann expansion are computed. They were aimed at improving the 
efficiency, accuracy or generality of the algorithms; the works of Yasui and Saikar^, and 
Fernandez Rico et are notable examples. Later, Maslen and Trefry 21 utilised an ap¬ 

proach based on the hypergeometric function which enabled to derive closed-form succinct 
analytical expressions for all necessary quantities. Despite those expressions were marred 
with numerical instabilities, it was a considerable progress at the time. Harris- 22 pursued 
the analytical approach of Maslen and Trefry, introduced considerable simplifications and 
several new expressions which allow more stable calculations of several auxiliary quantities. 

This paper is organised as follows. In Sec. [TT]we introduce the notation and recall relevant 
expressions from the previous works. In Sec. IIIII we introduce the large-order asymptotic 
expansion for the functions and verify the main results numerically. The corresponding 
expansion for the functions W tl is given in Sec. II VI Finally, in Sec. [V]we conclude our paper. 

II. PRELIMINARIES 

Let us consider a diatomic molecule placed on the z axis symmetrically around the origin. 
Slater-type orbitals (STOs) have the following generic form 

Xnim(r;C) = ^(Or^e-^^M), (1) 

where n and l are both integers such that n > l, (r, 9, 0) are the spherical coordinates of the 
given centre, S n (() = (2(0 n+1 / 2 / y/ (2 n)\ is the (radial) normalisation constant, and Y im are 
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spherical harmonics in the Condon-Shortley phase convention 

i i e im< ^ 

Yim( r) = QimPi {cos 0)—=, (2) 

v 27r 

where PJ 71 are the (unnormalised) associated Legendre polynomials^ and Vti m is the angular 
normalisation constant 




'2Z + 1(Z- M)! 


(l + \m\)\' 


( 3 ) 


Transformation to the real spherical harmonics, which are usually more convenient in cal¬ 
culations, can be performed with standard relations. 

Throughout the paper the electrons shall be denoted by 1,2,... and the nuclei by a, 6,.... 
All interparticle distances are shortly written as r pq , e.g., the distance between the first 
electron and the nucleus a is simply r\ a etc. (an exception is the internuclear distance for 
which the usual convention R := r ab is adapted). Let us introduce the prolate ellipsoidal 
coordinates, by means of the formulae 
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Via + r ib 

R 


Vi 


Via - r ib 

R 


( 4 ) 


where i = 1,2, and </>* are the corresponding azimuthal angles. The volume element becomes 
dvi = (|-) 3 (£ 2 — rff) d^i drji dfa. It is well known that the product of two Slater-type orbitals 
can be written in a closed-from in the prolate ellipsoidal coordinate system as follows 


!'i & 


K, 


Vi )Xn a l a m a ( r io,] Ca)Xn b l b m b ( r ib] C b) — 

r 


2tt 


■e-°*-** [(C 2 - 1)(1 - ril)] |Ml|/2 e iM >* S pq tf , 

p,q =o 


( 5 ) 


with M = m a — m b and T — l a + h + 2. The new coefficients are defined as a = — (£ 0 + ( b ), 


R\n a +n b +l 


The quantity 3 is a 


P = f («aCa + KbCfc), K Ob = S na {Ca)Sn b (Cb)^l a m a ^l b m b { f) 
square matrix with some numerical coefficients which can be tabulated. Details of this 

29] (see also the references therein). In conclusion, any 


transformation are given in Refs. 


nonzero two-centre electron repulsion integral over STOs can be written down as a finite 
linear combination of the following generic integrals 


f+i 


-2 n 


f+1 


-2tt 


!£%{*)= I / drn I dfa I d& 
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where explicit notation for the nonlinear parameters has been suppressed for brevity. The 
values of p t , qi and a are restricted to non-negative integers. 

Let us now introduce the Neumann expansion of the Coulomb interaction potential 


i o 00 ^ 

— = ±y y 

r\2 R ^ ^ 


-iy(2p + i) 


H) ! 

_(> + M) ! . 


0 <T =—(I 


(7) 


where £< = min(£i,£ 2 ) and £> = max(^,( 2 ), P™ are defined in the same way as in Eq. ([21), 
and Q™ are the associated Legendre functions of the second kind. By plugging the above 
expansion into Eq. ([ 6 ]) and after a straightforward integration over the angles one arrives at 

8 °° 

= x(- 1 ) a y( 2 V + 1 WZ(pi,P2,a 1 ,a 2 )il{qi,Pi)il(q 2 ,P2), ( 8 ) 

fl=<J 

where the basic quantities for the integration over 77 are 

= ( ~ 2 1)l ‘ j*' d,,pM(n)( 1 - pyip-e-f\ (9) 

and similarly for the £ integration 


W°{p l ,p 2 ,a 1 ,a 2 ) = w°(p 1 ,p 2 ,a 1 ,a 2 ) 

+ w^(p 2 ,pi,a 2 ,a 1 ), 


( 10 ) 


/ oo 

<% 1 q^i) y - lr^fe -^ 1 
x / <*6 p;(b) (& -1 

In general, the expansion given by Eq. (EJ) is infinite and terminates only in the special case 
of vanishing j3\ or /3 2 . Nonetheless, it is convergent for any physically acceptable values of 
the nonlinear parameters, i.e., oti > 0 and |/3j| < but the rate of convergence depends 
crucially on the values of fa. A practical observation is that larger values of /3; result in 
a slower convergence. Unfortunately, in actual calculations one can expect some of the 
integrals to approach the extreme case \[3,\ = cq. In such situation several tens of terms 
may be necessary to achieve convergence which significantly slows down the computations. 
Note parenthetically that the convergence is somewhat slower for larger values of a, but this 
effect is of secondary importance. 
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Calculation of the integrals i°(q,/3) is not connected with any significant overhead, even 
if large values of the parameters are necessary. Therefore, at present we see no reason to 
develop new methods for their computation. The available techniques appear to be entirely 
satisfactory and the recursive method put forward by Harris is particularly robust (see Ref. 


22] for an extended survey). We shall concentrate on the most difficult basic quantities, 


i.e., the integrals w°(pi,p 2 , ati, 012 )- Let us recall the analytical formula derived by Maslen 
and Trefry (after simplifications due to Harris) 


where 


w <7 Jpi,p 2 ,a 1 ,a 2 ) = 


(/1 + 0-)! 

2 - 

(p-a)l 



L' T (p 1 ,a 1 )kl(p2,a 2 ) 


P2+S 


s j =0 ^ ^ 


( 12 ) 


l Kp . a) = L-pLj f~ d( QZ(oe (( 2 - 1 r / v”«, (13) 

and 

k >- a > = |(rT^T I°° d( ~ < 14 ) 

The main goal of the present paper is to provide efficient and reliable methods for calculation 
of Wfj, for large values of /i. The problem can be solved in two ways. The first one is a direct 
attack by using the differential equation for W fl derived in the previous paper. The second 
method utilises Eq. (lT2]i and reduces the problem to calculation of L a (p, a) which appears to 
be more straightforward. In fact, large /j expansion of L°(p, a) is expected to be significantly 
less complicated than the corresponding one for W However, there is an additional cost 
of using Eq. (1121) which is absent in the first method where are calculated directly. Let 
us also note in passing that the auxiliary integrals k°(p,a), Eq. (1T4|) . can be computed 
efficiently with the available techniques and thus are not considered herein. 

Throughout the paper we rely on two special functions, E n (z) and a n (z). They are 
defined in the Appendix A and efficient methods of their computation are brieffy discussed. 
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III. LARGE-ORDER EXPANSION OF L^p- a) 


A. Initial reduction 


Let us recall two recursion relations which allow to simplify the problem significantly. 
They result directly from the properties of the Legendre functions and read 


L »(P+ !,«) 


(h + jjjjj+i(P, a) + oQ 

2 /i + 1 
2yU + 1 


(15) 

(16) 


By means of these recursions the necessary integrals L°(p-,a) can be efficiently computed 
starting with L^a) := Z°(0;a) only. Note that the above expressions require L^a) with 
even larger // than initially. In fact, they basically consist of increasing p and cr at cost of 
//. Therefore, the most important task is to calculate the integrals L^a) for large // with 
decent speed and precision. This is the main issue considered in the present section. 


B. Alternative integral representations of L^a) 


Our derivation starts with the differential equation for L /t (a) which was established in 


Ref. 


29| 


o?L n ^(oi) + 2aZ( i (cr) — \p(p + 1) + o 2 ] L^a) — —e a , (If) 

where the prime denotes differentiation with respect to a. Let us recall that the linearly 
independent solutions of the homogeneous differential equation are the well-known modified 
spherical Bessel functions^, i M (a) and /c M (a). This suggests that the desired solution of the 
inhomogeneous equation has the following form 

JC^a) i^a) +I M (a) k M (a), (18) 

where Z^a) and /C M (a:) are some functions which are yet to be determined. Let us addi¬ 
tionally enforce the constraint 


/C( t («) v(«) + Z^(a) k^a) = 0, 


(19) 
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valid strictly for every value of a > 0. Upon inserting the formula (lT8j) into the differential 
equation (ITT)) one obtains the following expression 


JC(a)z'(a)+l'(a)k'(a) = 


cr 


( 20 ) 


where we have taken advantage of the fact that i^a) and k^a) obey the homogeneous 
differential equation. The above expression and the constraint (1T9|) form the following system 
of linear equations 


( 21 ) 


Note that the determinant of the above 2x2 matrix (the Wronskian) is equal to — f ^2 
which is a direct consequence of the properties of the Bessel functions^. The system of 
linear equations (I2T1) can be solved right away, e.g., with the Cramer’s rule to give 


V(«) k^a) 




— e a /a 2 

^(«) k'^a) 


_w_ 


0 


K( a ) = + -e 

7r 

X 'M) = ~^ e_ %( a ), 


and after (indefinite) integration over a one arrives at 


( 22 ) 

(23) 


2 

7T 




dae a kJa) - kJa) 

7r 


dae °‘if i (a) 


(24) 


This expression is the general solution of the differential equation (ITT)) . In order to find 
a particular solution corresponding to the integrals (fT3l) we need to impose proper initial 
conditions. From Eq. (1TT1) one clearly sees t 


Additionally, the results presented in Ref. 


rat L^a) vanish as a —* oo for every /i. 
29| indicate that L^(a) vanish exponentially 
quickly in this limit [as e~° log (a) in the leading-order term]. This constitutes the first 
initial condition which we need to impose on the above general solution. The second initial 
condition results from the small a behaviour of L m (ck) 


= -7 e 






(2/x + l)!! (2/r + 1)!! 


log(2a) + A in + 0(a ), 


(25) 


where 7 e is the Euler-Mascheroni constant, and A in are some numerical coefficients inde¬ 
pendent of a (c./. the supplemental material to Ref. [291] ) - An essential feature of the above 
formula is the logarithmic singularity for small a which has to be reproduced by Eq. 
















The most succinct formula which takes both initial conditions into account reads 


L^a) = /C M (o) i M (a) + X f4 (a) k M (a), 


( 26 ) 


2 r 

M°0 = - 

^ J a 

2 r 

w = - / 

71 Jo 


/C M (a) = - dze Z k^(z), 


dze Z i^{z). 


(27) 


(28) 


Clearly, the formula (1261) is a new integral representation of the L^(a) functions, alternative 
to the definition given by Eq. (fT3jh 

At this point an extended comment is mandatory. One might be uncertain about the 
reason behind introduction of Eq. (1261) . It is clearly more complicated than the initial defi¬ 
nition, Eq. (fT3]h and appears to give no computational or theoretical advantages. However, 
it turns out that Eq. (TT3|) is a very inconvenient starting point for the present developments. 
Despite the large-order expansions of the Legendre functions, Q are well-known^ - —, they 
are too complicated to be used for our purposes. A naive approach where a large-order 
expansion of Q M is inserted into Eq. (TTH1) leads to intractable integrals requiring a numerical 
solution. On the other hand, Eq. (T26l) is formulated solely in terms of the modified spherical 
Bessel functions. This is advantageous, as the large-order expansions of i^a) and k^a) are 
more compact and straightforward. In particular, we rely on the recent works of Sidi and 
Hoggan^ 1 ^ where an elegant formulation has been given. For convenience of the readers we 
list the relevant formulae of Sidi and Hoggan in the Appendix B, utilising our notation. 

C. Large-order expansion of L^(a) 

Having the integral representation (|26|) at hand, it becomes straightforward to derive the 
large-order expansion of the pertinent integrals and /C^o). By inserting the integral 

representations (IBID and (1B2D into Eqs. (127D and (128D . respectively, one obtains 



m =0 


(30) 


(29) 
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after a simple exchange of variables. The coefficients in the expansions are obtained with 
help of Eq. (1B3|) 


/ I rn o 

dtEb m (at)e~ at = k ^( a /2) 2k a^ 2k (a), (31) 

k =1 

«W= f°° dt b -^-e~ M = jr,(-ir- k ^(a/2f t E^ 2t - 1 (a), (32) 

'' X k =1 

which is valid for m > 0. In the special case of m = 0 the corresponding results are 
Aq (a) = a M (a) and Aq(ck) = E^ + i(a). The modihed spherical Bessel functions in Eq. (EH 
which multiply the integrals X^(a) and /C M (cx) can also be expanded with help of Eqs. ED 
and (1B2[) . This leads to a product of two infinite series which can finally be rewritten as 


£ 


c s( a ) 


L " a 2^+1 + 0< + l/2)*’ 


k(a) + (-l)'Af (a) 


1=0 


(33) 

(34) 


A short remark on the mathematical nature of the above expansion is necessary. Note that 
Eq. (1331) would not be classified as an asymptotic expansion by some authors because the 
coefficients are explicitly /i-dependent. That is why we prefer to use the term large-order 
expansion. We verified that Eq. (I33|) is capable of providing arbitrarily accurate results if 
only the value of // is sufficiently large. Thus, from the pragmatic point of view, Eq. (|33ll 
gives an effective method to calculate the values of L^a) for large /i where other techniques 
run out of steam. 

From the point of view of some developments it is useful to analyse in details the first 
term of the expansion (j33j) . One easily arrives at 


LM = 


2/i + 1 


[£■„+! (a) + a„(a)] + h.o. 


(35) 


Additionally, if the large /j, asymptotic formulae for E^ + i(a) and a M (ct) are used, Eqs. 
and (IA7D . some simplifications occur and one finds 


L u (a) =- 

M ; 2/i + l 


1 1 
+ 


/i — a /i + a 




(36) 


provided that /i > a. 
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D. Numerical tests and examples 


It is now mandatory to verify how the new formula (1551) works in practice. In Table [Tj 
we present results of some exemplary calculations of L^a) with help of the new formula, 
Eq. (155|) . Different values of a and p are tested to find the actual range of applicability. 
Additionally, the number of terms in the infinite expansion (1551) necessary to reach the 
maximal possible precision was listed in each case. A more detailed inspection of Table Q] 
reveals some general conclusions about the range of the parameters where Eq. (155|) gives 
sufficiently accurate results. One sees that the convergence of the infinite summation in 
Eq. (1551) is excellent for small or moderate values of a. Unfortunately, it deteriorates 
quickly when the values of a and p approach each other. In the case when a > p no useful 
information about L^(a) can be obtained with help of Eq. (155]) . However, this is not a 
reason for a major concern. In fact, the large a expansion of L^(a) was given in Ref. 29] 
and it works reasonably well for both small and large values of p. We conclude that Eq. 
(1551) is a preferred computational technique when p is large and a is small or moderate at 
the same time. 


IV. LARGE-ORDER EXPANSION OF W^{ Pl , p 2 ; a 1} a 2 ) 
A. Initial reduction 


Let us reduce the number of independent parameters in the integrals W fl by using two 
convenient formulae. The first one is the remainder in the recursive method proposed by 
Kotani^ 


^ +1 (Pi,P2,ai,a 2 ) = —— ~~~~y + ^ W£+i(Pi,P 2, Oh, <* 2 ) - {p - a)(p + a + 1) 

x w;(p! + i,p 2 + 1 , a 2 ) + ^ + a ^ 1 l^ + a)2 K-M’P*i « 2 ). 


(37) 


This expression is numerically stable for a wide range of the parameters values. As a result, 
it constitutes a reliable method for computation of W°(pi,p 2 , «i, a 2 ) from the integrals with 
a = 0. Additionally, the values of p 2 can be increased by differentiation 


W°(pi,p 2 ,a 1 ,a 2 ) 


f)P2 

' 1)B a5 5jW " (pi ’0’“ 1 ’“ 2 )' 


(38) 
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Further in the article, we consider the large /i expansion of the basic integrals W p (p \, 0, op, 0 : 2 ). 
Note that differentiation with respect to an could be used to increase the value of pi, but 
this approach is not particularly advantageous in the present context. For convenience, we 
introduce the following shorthand notation, ai, a 2 ) = PF°(p, 0 , ai, a 2 ). 


B. Alternative integral representations of W p (pi, 0, ai, a 2 ) 


Let us recall the differential equation for W p (p\ an, a 2 ) obtained in Ref. 



0^ d 

a 2 ai, a 2 ) + 2a 2 —— W^p; ai, a 2 ) + 


/i(/i + 1) + a^ ai,a 2 ) — —E_ p (ai + a 2 ), 


(39) 


which provides the starting point for our derivation. Note that the solutions of the homo¬ 
geneous equation are well-known and are the same as for Eq. cm. Therefore, the solution 
can be written in the form analogous to Eq. (TT 8 T) and the derivation follows along a very 
similar line as for L p . There is no need to repeat details of the derivation and we present 
only the final result 


ai, a 2 ) — /C M (p; aq, a 2 ) i M (a 2 ) + X M (p; ai, a 2 ) fc A 1 (a 2 ), 

2 

/C Al (p; ai, a 2 ) = - / dz k p {z) E_ p (a x + z), ( 40 ) 

2 f a2 

l p (p;a 1 ,a 2 ) = - / dzi^z) E_ p (a 1 + z), 

^ Jo 


by imposing proper initial conditions (c./. Ref. 29]). Note that the basic integrals were 
expressed though the modified spherical Bessel functions, in analogy with L p functions 
considered before. Clearly, Eqs. (FIU1) may be useful on their own ( e.g. evaluation by 
a numerical integration), but in the present paper we concentrate solely on the large p 
expansion of W p {jp ai, a 2 ). 
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C. Large-order expansion of W°(pi, 0, «i, 02 ) 


Let us now insert the asymptotic expansions of i M (ct) and kp(a), Eqs. (IBID and (IB2D . 
into the integral representation (l40lh After straightforward rearrangements one arrives at 


2 /a 2 \ m +! 


i> ;ai , Q2 ) = -=(^) 


T,T(ai.a2) 
r(M + 3 / 2 )^ 0 ' + l/ 2 )’“’ 


E 


1 itM, 

h (/< + i/2) m 


/C M (p;ai,a 2 ) = —j= ( — ) T(/i + l/2)V 


-i) m , 


(41) 

(42) 


which is an analogue of Eqs. (l29]h The analytical formulae for the coefficients are obtained 
by recalling Eq. 


r^f(ai,a 2 )= / dtt^ b m (a 2 t) E_ p (ai + a 2 t) 


EM) 


k =1 


m-k Sink 

k\ 


(43) 


(a 2 /2) 2k u M + 2 fe,p( a; i) a 2 ), 


TT(a i,a 2 ) 


J dt b E_ p (ai + a 2 t ) 
E ( _ ir -*ir (“2/2)“ c 2 ), 


(44) 


for m > 1, and T^(a \, ct 2 ) = a7 /ip (ai, ct 2 ), T^(a\, a 2 ) = fi_ M _i )P (ai, ct 2 ). The basic integrals 
are defined as 


(an,ay) = / dtt n E_ p (ai + a 2 t), 


f2rap(«i, a 2 ) = / dtt n E_ p (ai + a 2 f). 


(45) 

(46) 


Note that evaluation of Eqs. (1431) and (1441) requires (n np with n > 0, p > 0, but in the 
case of Q np the values of n can be negative. Additionally, the Erst argument (n) in u> np is 
always larger than /i [c.f. Eq. (145D ], The present method is intended to be used for large 
/i and we concentrate on evaluation of u np with large n. Unfortunately, for the integrals 
Q np such simplihcations do not occur and more general methods are required. Calculation 
of the basic integrals is discussed in the next section, with a considerable emphasis on the 
numerical stability. 
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Finally, one combines the asymptotic expansions (jUJ) with the initial formula, Eq. 
and after some rearrangements the following expression is obtained 

^ OO 


W M (p; ai,a 2 ) = 


£ 


d" 


2^ + i^(^ + i/2r’ 


where we have suppressed the notation for the nonlinear parameters, and 

S 

df = Y (— I i''' y‘(a 1 ,a 2 ) + (-l)*rf p (ai,Q:2) 


/=o 


( 47 ) 


(48) 


which constitutes the main result of the present section. 


D. Calculation of the basic integrals 

Let us begin with calculation of the integrals oj np . Integration of Eq. (|45|) by parts leads 
to the following recursion 

1 (^2 

^np(aq, 012 ) = —— E_ p (a\ + 02 ) H- — u n + i,p+i(o!i + ot2t). (49) 

n + 1 n + 1 

In principle, the above relation can be used to calculate the values of oo np by downward 
recursion, starting at some large n with an arbitrary value. However, the main drawbacks 
of this approach are difficulties in controlling the error and choice of the starting point. 
Therefore, we propose to iterate this recursion analytically N times which gives 

N u 

^n P {oii, CK 2 ) = n\ y' 7 — , 2 1 \, E_ p _ k {ai + ck 2 ) + Rn, (50) 

('n + k + l)\ 

Rn = w n+ jv+i, P +jv+i(o!i, ot2). (51) 

Note that the above expression is formally exact for each N. Additionally, when n is large 
the terms in the above sum vanish very quickly and large values of N give very small 
contributions to the total value of the integral. Similarly, the remainder Rn vanishes fast 
with increasing N. To estimate in advance the required values of N we establish approximate 
upper bounds for the values of R n (note that Rn is positive by definition). Let us first 
insert the integral representation of E_ p , Eq. flAljl . into Eq. (145| and reverse the order of 
integrations. One arrives at the alternative integral representation of the remainder 

/ OO /* 1 

dzz p+N+1 e~ aiz J dtt n+N+1 e~ a2tz , (52) 
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which is strictly bounded from above by 

POO p 1 -| 

Rn< dzz p+N+1 e- a ' z dtt n+N+1 = ———E_ p _ N _ 1 (a 1 ). (53) 

J 1 Jo 71 + N + 2 

Additionally, one can verify that £L p _jv-i(ai) is bounded from above by (p+ N +1)!/ a p+N+1 . 
This finally gives the estimation 


„ n\ ( a 2 V v+1 (p + N + 1)! 

' (n + JV + 2)f 

Passing to the integrals fi np , the optimal algorithm depends on the sign of n. 
as before, by inserting the integral representation of E_ p , Eq. (lAlll . into Eq. 
exchanging the order of integrations 


(54) 

Similarly 
jSD and 


H np (a; i,a 2 )= / dzz p e 


— Ot\Z 


dtt n e 


n —OL 2 tz _ 


dzz p e~ aiz E- n (a 2 z). (55) 


When the values of n are positive one can use the recursive relation (1A3I) which gives 

72 1 

^np(o ; l, CX2) — - ^n-l,p-l(ttl, CK2) H - E-p+l( a l + a 2)- 

a 2 a 2 


(56) 


This recursion relation is completely stable when carried out in the upward direction along 
the “diagonal” lines. However, it is not self-starting and requires values of Q n0 and Hop to 
initiate. Analytical expression for the latter is fairly obvious, H 0p = E_ p+ i(ai + a 2 ), and 
calculation of the former is based on the following relation 

C~ ai CM, 1 

H n o( a l; a 2) = - E-n+l{o! 2 ) -H n _ li o(o:i, Ot 2 )i (57) 

a 2 oi 2 


starting with H 00 = Ei(a.\ + a 2 ). This recursion is stable provided that the value of a 2 is 
moderate or large. If a 2 is small the following series expansion is used 


0 / n _ n '- p ( n v- (~a 2 ) k E_ k (ati) 

H n o( tt l5 a 2) — n+ i En+ l(oq) — 


a 


^ k\ n + k + 1 ’ 

fc =0 


(58) 


which can be derived by using elementary methods. Similarly as before one can verify that 
the absolute value of each term in the above sum is bounded by a^/oJ^ 1 which can be used 
to estimate the convergence rate. 

Finally, let us discuss calculation of Vt np for negative values of n. The following recursion 
can be derived with help of Eqs. (I46|) and 


H np (o( 1, O 2 ) 


1 , 0 ( 2 ) ——H n _ 1)P (o(i, a 2 ) H- E\_ n (cx 2 ), 

Oi 2 a 2 a 2 


(59) 


15 












which can be used to increase p at cost of n. This recursive relation introduces some 
instabilities into the calculation, but this fact is not significant as the values of p rarely 
exceed 10. To initiate the above recursion one requires the values of fi n0 . Similarly as 
before, the following expression is straightforward to derive 

6 _ai CX\ 

^no( a l; a 2 ) = - Ei _ n (ay)-^n—l,o( a l) Qg). (60) 

Ot2 O' 2 

This recursion relation is carried out downward, starting with tt-N,o at some large N. This 
completes the formalism of calculation of the basic integrals. 


E. Numerical tests and examples 


In Tables El and IIIII we present results of exemplary calculations of the W^{p] ay, ay) 
functions with Eq. (1471) and comparison with the exact results. One can see that the 
method based on Eq. (147|) converges in at most few tens of terms provided that ay and ay 
are both small or moderate. In fact, for small values of the nonlinear parameters we managed 
to obtain the convergence even for p as small as 10 which shows the potential of the method. 
Unfortunately, when ay and a 2 are both large (larger than 50, say) the series (TTTT) have an 
oscillatory behaviour and no convergence was achieved after summing 200 terms. However, 
in the regime of large ay and ay one can resort to different techniques e.g. asymptotic 


expansions presented in Ref. 


29]. Moreover, when ay and ay are simultaneously large the 


resulting integrals are expected to be very small and they are likely to be negligible. Table 
nn lists the results for p = 0 whilst the corresponding values for p = 8 are given in Table 
EH A more detailed comparison reveals that larger values of p are connected with slower 
convergence of the series fjTTjl . but the range of applicability remains roughly the same. 


V. CONCLUSIONS 

We have presented a new systematic approach to the calculation of basic quantities ap¬ 
pearing in the ellipsoidal expansion of the two-electron integrals over Slater-type orbitals. 
Large-order (//) expansions of the functions L /t and have been given and their accuracy 
and range of applicability has been determined numerically. The new method allows to 
calculate higher-order terms of the Neumann expansion with a significantly reduced com¬ 
putational cost. As a result, this is a step towards reduction of the gap in computational 


16 







timings between STOs and GTOs. Moreover, the presented expressions may be useful in 
mathematical studies of convergence of the Neumann series and rational design of conver¬ 
gence acceleration techniques. 
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Appendix A: Auxiliary integrals 


Virtually all final working formulae obtained in the present paper are given in terms of 
the basic integrals E n (z) and a n (z). They are defined through the integral representations 

/ oo -zt 

dt—, (Al) 

a n (z) = [ dtt n e~ z \ (A2) 

Jo 

where n is an arbitrary integer in the former and a nonnegative integer in the latter. Cal¬ 
culation of a n is most easily carried out with help of the Miller algorithm 49 - as discussed by 
Harris 22 . The integral E n is usually called the generalised exponential integral. Computa¬ 
tion of E n differs depending on the sign of n. For a negative integer n the following recursion 
is completely stable in the upward direction 

E n (z) = — — E n _i(z) + —. (A3) 

z z 


For positive n and z < 1 one uses the series expansion 


-z 


\TL—1 


En ^ Z ’ (n — 1)! 


[4/ (n) — log z] — 


-zY 


k=0 
k^n—1 


k\ (1 — n + k) ’ 


(A4) 


where 'k(n) is the digamma function at integer argument. The above infinite summations 
converge to the machine precision in, at most, few tens of terms. Finally, for positive n and 
z > 1 the continued fraction (CF) formula can be applied 


E n {z) = e z 


/lplp+12 \ 

Vi+1+I+ 1+ z+'"J ' 


(AS) 
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To evaluate the CF one can use the Lentz algorithm^. The only inconvenience is that 
consecutive numerators and denominators in the Lentz scheme grow very quickly with the 
number of terms retained in Eq. (I A 511 . Therefore, it is necessary to rescale them from time 
to time by a small number to avoid numerical overflows. Let us also recall the leading terms 
of the large-order asymptotic expansions for E n {z) and a n (z) which read 

E n (z) = -^- + o(-\, (A6) 

n + z \ n J 

aJz) = — + o(E], (A7) 

71 — Z \ n J 

where n > z. 


Appendix B: Large-order asymptotic formulae for i^a) and k^{a) 


According to the work of Sidi et a/. 47 ' 48 the modified spherical Bessel functions posses the 
following large-order expansions 


\A" (*/2) M b m (z) 

— r 0“+t)ii> + i/ 2 ) m ’ 

_ r(h + j) \ ^ b m (z) , - 

4 (z/2)^ +1 (ji + l/2) m 


(Bl) 

(B2) 


as /i —» 00 at a hxed z, where 


m 


k= 1 


(B3) 


for m > 0 and 5 0 (a) = 1- The quantities S m k in the above expression are the Stirling 
numbers of the second kind4^ defined recursively as 


SmO UmCh S m l 1, S mm 1, 

*S 'mk l,fc—1 T kS m —i 


and we additionally adapt the convention S m k = 0 for m < k or m < 0. 


(B4) 

(B5) 
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TABLE I: The functions L^a) calculated for some representative values of a and //. Exact denotes 


values calculated using explicit expressions (c./. Ref. 


291 ]) in extended arithmetic precision with 


the Mathematica package (all digits shown are correct). Large-order expansion column shows 
results of calculations with Eq. (|47fl in the double precision arithmetic. Convergence denotes a 
number of terms in Eq. (1471) required to converge the summation to relative precision of 2 • 10~ 16 . 
The symbol [k] denotes the powers of 10, 10 fc . 



exact 


large-order expansion 

convergence 

a = 0.1 

30 

9.72 733 864 877 071 

[-04] 

9.72 733 864 877 071 

[-04] 

9 

40 

5.51 662 783 117 224 

[-04] 

5.51 662 783 117 224 

[-04] 

9 

50 

3.54 810 355 237 372 

[-04] 

3.54 810 355 237 372 

[-04] 

7 

60 

2.47 209 822 882 328 

[-04] 

2.47 209 822 882 328 

[-04] 

7 

a = 1.0 

30 

3.94 720 438 208 518 

[-04] 

3.94 720 438 208 518 

[-04] 

13 

40 

2.24 043 509 438 319 

[-04] 

2.24 043 509 438 319 

[-04] 

11 

50 

1.44 153 386 177 520 

[-04] 

1.44 153 386 177 520 

[-04] 

9 

60 

1.00 458 613 132 488 

[-04] 

1.00 458 613 132 488 

[-04] 

9 

a = 10.0 

30 

4.78 078 398 572 794 

[-08] 

4.78 078 398 572 794 

[-08] 

29 

40 

2.73 528 592 642 051 

[-08] 

2.73 528 592 642 051 

[-08] 

21 

50 

1.76 662 929 639 839 

[-08] 

1.76 662 929 639 839 

[-08] 

19 

60 

1.23 372 624 613 915 

[-08] 

1.23 372 624 613 915 

[-08] 

19 

a = 30.0 

30 

9.48 264 212 820 654 

[-17] 

divergence 


— 

40 

5.51 072 668 384 366 

[-17] 

divergence 


— 

50 

3.58 705 498 786 413 

[-17] 

3.58 705 498 786 412 

[-17] 

51 

60 

2.51 610 451 657 133 

[-17] 

2.51 610 451 657 133 

[-17] 

41 


21 















TABLE II: The functions ILL (p: aq, 0 : 2 ) calculated for some representative values of the parameters. 


Exact denotes values calculated using explicit expressions ( c.f. Ref. 


29l |) in extended arithmetic 


precision with the Mathematica package (all digits shown are correct). Large-order expansion 
column shows results of calculations with Eq. (1331) in the double precision arithmetic. Convergence 
denotes a number of terms in Eq. (1331 required to converge the summation to relative precision of 
2 • 1CT 16 . The symbol [k] denotes the powers of 10, 10 fc . 



exact 


large-order expansion 

convergence 

aq = 1.0, 0:2 = 1.0, p = 0 

30 

7.26 438 420 525 738 

[-05] 

7.26 438 420 525 741 

[-05] 

13 

40 

4.12 230 721 703 987 

[-05] 

4.12 230 721 703 989 

[-05] 

11 

50 

2.65 207 347 050 162 

[-05] 

2.65 207 347 050 163 

[-05] 

9 

60 

1.84 808 542 738 505 

[-05] 

1.84 808 542 738 506 

[-05] 

7 

aq = 1.0, 012 = 5.0, p = 0 

30 

4.43 294 559 944 128 

[-07] 

4.43 294 559 944 129 

[-07] 

19 

40 

2.51 607 499 568 558 

[-07] 

2.51 607 499 568 559 

[-07] 

16 

50 

1.61 886 531 395 923 

[-07] 

1.61 886 531 395 924 

[-07] 

12 

60 

1.12 815 856 115 223 

[-07] 

1.12 815 856 115 224 

[-07] 

9 

ctq = 1.0, 0:2 = 10.0, p = 0 

30 

1.62 914 653 440 044 

[-10] 

1.62 914 653 440 043 

[-10] 

29 

40 

9.24 696 775 490 950 

[-10] 

9.24 696 775 490 946 

[-10] 

21 

50 

5.94 963 334 102 345 

[-10] 

5.94 963 334 102 343 

[-10] 

19 

60 

4.14 621 345 113 117 

[-10] 

4.14 621 345 113 115 

[-10] 

19 

oq = 1.0, 0:2 = 50.0, p = 0 

30 

1.49 278 026 722 289 

[-27] 

divergence 


- 

40 

8.47 302 548 103 673 

[-28] 

divergence 


- 

50 

5.45 168 822 712 533 

[-28] 

divergence 


- 

60 

3.79 921 108 103 191 

[-28] 

3.79 921 108 103 191 

[-28] 

157 

aq = 5.0, aq = 5.0, p = 0 

30 

4.85 315 094 815 931 

[-09] 

4.85 315 094 815 933 

[-09] 

21 


Continued on the next page 
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TABLE II: Continued from the previous page. 


40 

2.75 

906 

112 

497 

927 [- 

-09] 

2.75 

906 112 497 928 [ 

-09] 

17 

50 

1.77 

656 

445 

608 

724 [- 

-09] 

1.77 

656 445 608 724 [ 

-09] 

15 

60 

1.23 

857 

631 

621 

439 [- 

-09] 

1.23 

857 631 621 440 [ 

-09] 

13 






Ot\ 

= 5.0 

'> a 2 = 

10.0, p = 0 



30 

2.17 

631 

353 

217 

224 [- 

-11] 

2.17 

631 353 217 225 [ 

-11] 

29 

40 

1.23 

815 

626 

851 

864 [- 

-11] 

1.23 

815 626 851 865 [ 

-11] 

21 

50 

7.97 

525 

834 

542 

020 [- 

-12] 

7.97 

525 834 542 023 [ 

-12] 

19 

60 

5.56 

120 

404 

111 

395 [- 

-12] 

5.56 

120 404 111 398 [ 

-12] 

19 






Ot\ 

= 5.0 

, OL2 = 

50.0, p = 0 



30 

2.51 

589 

622 

982 

829 [- 

-29] 


divergence 


- 

40 

1.43 

273 

896 

436 

017 [- 

-29] 


divergence 


- 

50 

9.23 

286 

988 

543 

448 [- 

-30] 


divergence 


- 

60 

6.43 

979 

747 

051 

030 [- 

-30] 

6.43 

979 747 051 030 [ 

-30] 

189 






Ot\ 

— 10.0, 012 — 

10.0, p = 0 



30 

1.09 

589 

723 

489 

100 [- 

-13] 

1.09 

589 723 489 100 [ 

-13] 

29 

40 

6.24 

425 

513 

575 

906 [- 

-14] 

6.24 

425 513 575 909 [ 

-14] 

21 

50 

4.02 

496 

564 

881 

221 [- 

-14] 

4.02 

496 564 881 222 [ 

-17] 

19 

60 

2.80 

774 

946 

310 

143 [- 

-14] 

2.80 

774 946 310 144 [ 

-17] 

19 






Ot\ 

— 10.0, 012 — 

50.0, p = 0 



30 

1.54 

163 

263 

060 

580 [- 

-31] 


divergence 


- 

40 

8.80 

879 

102 

228 

926 [- 

-32] 


divergence 


- 

50 

5.68 

572 

697 

434 

438 [- 

-32] 


divergence 


- 

60 

3.96 

925 

492 

004 

272 [- 

-32] 


divergence 


- 






Ot\ 

= 50.0, oi2 = 

50.0, p = 0 



30 

3.80 

314 

546 

909 

033 [- 

-49] 


divergence 


- 

40 

2.20 

243 

367 

882 

466 [- 

-49] 


divergence 


- 

50 

1.43 

106 

434 

831 

316 [- 

49] 


divergence 


- 

60 

1.00 

278 

875 

314 

209 [- 

-49] 


divergence 


- 


23 















TABLE III: The functions WRp; aq, a.^) calculated for some representative values of the param¬ 


eters. Exact denotes values calculated using explicit expressions ( c.f . Ref. 


291 ]) in extended 


arithmetic precision with the Mathematica package (all digits shown are correct). Large-order 
expansion column shows results of calculations with Eq. (1331) in the double precision arithmetic. 
Convergence denotes a number of terms in Eq. (1331) required to converge the summation to relative 
precision of 2 • lCP 16 . The symbol [k] denotes the powers of 10, 10 fc . 



exact 


large-order expansion 

convergence 

aq = 1.0, 0:2 = 1.0, p = 8 

30 

8.57 988 797 367 550 

[-02] 

8.57 988 797 367 554 

[-02] 

13 

40 

4.83 735 222 300 728 

[-02] 

4.83 735 222 300 730 

[-02] 

11 

50 

3.10 265 767 847 292 

[-02] 

3.10 265 767 847 291 

[-02] 

9 

60 

2.15 848 287 860 469 

[-02] 

2.15 848 287 860 470 

[-02] 

9 

aq = 1.0, 0.2 = 5.0, p = 8 

30 

3.76 412 026 039 309 

[-06] 

3.76 412 026 039 311 

[-06] 

19 

40 

2.10 468 765 565 305 

[-06] 

2.10 468 765 565 306 

[-06] 

17 

50 

1.34 482 396 704 125 

[-06] 

1.34 482 396 704 126 

[-06] 

15 

60 

9.33 659 812 480 290 

[-07] 

9.33 659 812 480 295 

[-07] 

13 

ctq = 1.0, «2 = 10.0, p = 8 

30 

4.37 682 376 980 937 

[-09] 

4.37 682 376 980 939 

[-09] 

29 

40 

2.45 407 770 819 367 

[-09] 

2.45 407 770 819 369 

[-09] 

21 

50 

1.57 011 062 044 134 

[-09] 

1.57 011 062 044 135 

[-09] 

19 

60 

1.09 084 224 914 003 

[-09] 

1.09 084 224 914 004 

[-09] 

19 

oq = 1.0, ot .2 = 50.0, p = 8 

30 

1.79 779 640 123 314 

[-27] 

divergence 


— 

40 

1.01 180 529 453 933 

[-27] 

divergence 


— 

50 

6.48 466 217 722 202 

[-28] 

divergence 


- 

60 

4.50 946 436 202 621 

[-28] 

4.50 946 436 202 621 

[-28] 

181 

aq = 5.0, aq = 5.0, p = 8 

30 

1.44 602 672 357 563 

[-08] 

1.44 602 672 357 562 

[-08] 

21 


Continued on the next page 
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TABLE III: Continued from the previous page. 


40 

8.19 

252 

305 

004 

658 [- 

-09] 

8.19 

252 305 004 661 [ 

-09] 

17 

50 

5.26 

663 

146 

433 

895 [- 

-09] 

5.26 

663 146 433 898 [ 

-09] 

15 

60 

3.66 

849 

871 

467 

375 [- 

-09] 

3.66 

849 871 467 377 [ 

-09] 

13 






Ol\ 

= 5.0 

1 «2 = 

10.0, p = 8 



30 

4.24 

176 

585 

818 

655 [- 

-11] 

4.24 

176 585 818 656 [ 

-11] 

29 

40 

2.40 

082 

417 

196 

577 [- 

-11] 

2.40 

082 417 196 578 [ 

-11] 

21 

50 

1.54 

267 

805 

985 

739 [- 

-11] 

1.54 

267 805 985 739 [ 

-11] 

19 

60 

1.07 

428 

909 

373 

431 [- 

-11] 

1.07 

428 909 373 431 [ 

-11] 

19 






Ol\ 

= 5.0 

1 «2 = 

50.0, p = 8 



30 

2.98 

140 

021 

972 

450 [- 

-29] 


divergence 


- 

40 

1.68 

614 

952 

443 

023 [- 

-29] 


divergence 


- 

50 

1.08 

306 

875 

847 

441 [- 

-29] 


divergence 


- 

60 

7.54 

082 

251 

283 

390 [- 

-30] 


divergence 


— 






Ot\ 

— 10.0, 012 — 

10.0, p = 8 



30 

1.75 

800 

825 

092 

494 [- 

-13] 

1.75 

800 825 092 495 [ 

-13] 

29 

40 

9.99 

225 

776 

925 

089 [- 

-14] 

9.99 

225 776 925 094 [ 

-14] 

21 

50 

6.43 

335 

525 

002 

640 [- 

-14] 

6.43 

335 525 002 643 [ 

-14] 

19 

60 

4.48 

491 

212 

960 

371 [- 

-14] 

4.48 

491 212 960 374 [ 

-14] 

19 






Ot\ 

— 10.0, 012 — 

50.0, p = 8 



30 

1.79 

687 

571 

912 

046 [- 

-31] 


divergence 


- 

40 

1.02 

no 

421 

809 

320 [- 

-31] 


divergence 


- 

50 

6.57 

355 

034 

007 

000 [- 

-32] 


divergence 


- 

60 

4.58 

239 

720 

894 

004 [- 

-32] 


divergence 


- 






Ot\ 

= 50.0, oi2 = 

50.0, p = 8 



30 

4.14 

652 

002 

770 

361 [- 

-49] 


divergence 


— 

40 

2.39 

745 

614 

502 

924 [- 

-49] 


divergence 


— 

50 

1.55 

652 

302 

599 

839 [- 

49] 


divergence 


— 

60 

1.09 

019 

745 

763 

500 [- 

-49] 


divergence 


- 


25 
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